2 Codigo Reproducible: Imputación Múltiple - Recalibración de modelo KFRE para predecir falla renal en asegurados de EsSalud

1 Code to reproduce results of the manuscript ‘Kidney Failure Prediction: Multicenter External Validation of the KFRE Model in Patients with CKD Stages 3-4 in Peru’

1.1 Introduction

This document presents the code and results of the main analysis shown in the article.

1.2 Setup

rm(list = ls())

# Use pacman to check whether packages are installed, if not load
if (!require("pacman")) install.packages("pacman")

Loading required package: pacman

library(pacman)

# Unload all package to begin in a session with only base packages
pacman::p_unload("all")

The following packages have been unloaded:
pacman

# Install packages
pacman::p_load(
  here, 
  skimr, 
  survival,
  rms,
  cmprsk,
  riskRegression,
  mstate,
  pseudo,
  pec,
  riskRegression,
  plotrix,
  knitr,
  splines,
  kableExtra,
  flextable,
  gtsummary,
  boot,
  tidyverse,
  rsample,
  gridExtra,
  webshot, 
  patchwork,
  survival, 
  cmprsk, 
  survminer, 
  ggsci, 
  cowplot, 
  scales, 
  patchwork, 
  labelled, 
  glue, 
  dcurves, 
  broom, 
  downlit, 
  xml2, 
  gghalves, 
  devtools, 
  htmltools, 
  gghalves, 
  ggtext, 
  DiagrammeR, 
  gt, 
  janitor, 
  VIM, 
  PerformanceAnalytics, 
  mice, 
  rms, 
  naniar, 
  DescTools, 
  gtools, 
  ggExtra, 
  furrr, 
  future, 
  tictoc, 
  parallel,
  ggmice
)

if (!require("impstep")) remotes::install_github("bgravesteijn/impstep", force = TRUE)

Loading required package: impstep

if (!require("smplot2")) devtools::install_github('smin95/smplot2', force = TRUE)

Loading required package: smplot2

Updated tutorial for smplot: smin95.github.io/dataviz/

# You will need Rtools to install packages from Github on Windows
# `devtools` with throw an informative error if Rtools is not found
# if (!"devtools" %in% installed.packages()) install.packages("devtools")
# devtools::install_github("jesse-smith/futuremice")

library(impstep)

## Revisar:; https://amices.org/ggmice/index.html

1.3 Cargar datos

Los datos completos se muestran a continuación, luego de seleccionar a las variables con las que trabajaremos:

# Import data
data <- readRDS(here::here("Data", "Tidy", "datos_total_integrados.rds")) |> 
  select(cas, sex, age, hta, dm, crea, 
         ckd_stage, ckd_stage2, 
         eGFR_ckdepi, acr, urine_album, 
         urine_crea, time5y, eventd5y, 
         grf_cat, acr_cat, ckd_class,
         death2y, eventd2ylab, death5y, 
         eventd5ylab, eventd, time) 

data |> 
  glimpse()

Rows: 137,475
Columns: 23
$ cas         <chr> "Arequipa", "Arequipa", "Arequipa", "Arequipa", "Arequipa"…
$ sex         <fct> Masculino, Masculino, Masculino, Masculino, Masculino, Fem…
$ age         <dbl> 15, 15, 15, 15, 15, 18, 15, 16, 21, 17, 14, 15, 14, 14, 15…
$ hta         <dbl> 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1…
$ dm          <dbl> 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0…
$ crea        <dbl> 0.99, 1.20, 0.97, 0.91, 0.81, 0.76, 0.88, 1.20, 0.61, 1.07…
$ ckd_stage   <fct> Stages 1-2 y 5, Stages 1-2 y 5, Stages 1-2 y 5, Stages 1-2…
$ ckd_stage2  <fct> Stages 1-3 y 5, Stages 1-3 y 5, Stages 1-3 y 5, Stages 1-3…
$ eGFR_ckdepi <dbl> 113.08736, 89.62041, 115.91243, 125.21489, 132.51473, 114.…
$ acr         <dbl> NA, NA, NA, NA, NA, NA, 0.2522603, NA, 392.8585677, NA, NA…
$ urine_album <dbl> NA, NA, NA, NA, NA, NA, 22.60, NA, 196.94, NA, NA, NA, NA,…
$ urine_crea  <dbl> NA, NA, NA, NA, NA, NA, 89.5900, NA, 0.5013, NA, NA, NA, N…
$ time5y      <dbl> 5.0000000, 5.0000000, 5.0000000, 5.0000000, 5.0000000, 4.0…
$ eventd5y    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ grf_cat     <fct> G1, G2, G1, G1, G1, G1, G1, G2, G1, G1, G1, G1, G2, G2, G1…
$ acr_cat     <fct> NA, NA, NA, NA, NA, NA, A1, NA, A3, NA, NA, NA, NA, NA, NA…
$ ckd_class   <fct> NA, NA, NA, NA, NA, NA, Low risk, NA, High risk, NA, NA, N…
$ death2y     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ eventd2ylab <chr> "Alive w/o Kidney Failure", "Alive w/o Kidney Failure", "A…
$ death5y     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ eventd5ylab <chr> "Alive w/o Kidney Failure", "Alive w/o Kidney Failure", "A…
$ eventd      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ time        <dbl> 7.4934976, 7.6933607, 7.4934976, 7.9041752, 7.9041752, 4.0…

1.4 Análisis inicial de datos

Vamos a identificar datos de tiempo no plausibles. Se aprecia que todos los datos de tiempo que no son plausibles pertenecen a individuos que hicieron falla renal, y esta implausibilidad se debe a tiempos hasta eventos negativos. Respecto a los datos perdidos, vemos que hay datos perdidos de tiempo tanto en individuos vivos que no reportaron falla renal como en individuos que reportaron falla renal. Como era de esperarse, no hubo datos perdidos de tiempo en individuos que fallecieron, dado que estos datos provienen de la oficina de asegurados que cruza datos con RENIEC. Asimismo, un grupo importante de individuos tuvo datos perdidos de tiempo y de evento (no se sabe si desarrollaron o no falla renal o muerte).

data |> 
  mutate(eventdf = as.character(eventd),
         eventdf = case_match(eventdf, 
                              "0" ~ "Vivo y sin falla renal", 
                              "1" ~ "Falla renal", 
                              "2" ~ "Muerto sin falla renal", 
                              NA ~ "Dato perdido")) |> 
  bind_rows(data |> mutate(eventdf = as.character("Total"))) |> 
  mutate(eventdf = factor(eventdf, 
                          levels = c("Vivo y sin falla renal", 
                                     "Falla renal", 
                                     "Muerto sin falla renal", 
                                     "Total",
                                     "Dato perdido"))) |> 
  count(time, eventdf) |> 
  mutate(time = case_when(is.na(time) ~ 15, 
                          TRUE ~ time), 
         time_tipo = case_when(time > 0 & time < 15 ~ "Valor Plausible", 
                               time <= 0 ~ "Valor Implausible", 
                               time == 15 ~ "Dato perdido", 
                               TRUE ~ as.character(NA))) |> 
  ggplot(aes(x = time, y = n, color = time_tipo)) +
  geom_point(shape = 21, alpha = 0.5) + 
  geom_vline(xintercept = 0.1, linetype = 2, color = "red") +
  scale_y_continuous(trans = "log10") + 
  labs(color = "Tipo de dato") + 
  facet_wrap(. ~ eventdf) + 
  theme_bw() -> p_tiempo_perdido

ggsave(filename = "p_tiempo_perdido.png", 
       plot = p_tiempo_perdido, 
       device = "png", 
       path = here("Figures", "Imputation_Diagnostic"), 
       scale = 2.5, 
       width = 9,
       height = 9, 
       units = "cm", 
       dpi = 600, 
       bg = "white") 

Respecto a la distribución de estos datos perdidos, tenemos lo siguiente:

data |> 
  mutate(eventdf = as.character(eventd),
         eventdf = case_match(eventdf, 
                              "0" ~ "Vivo y sin falla renal", 
                              "1" ~ "Falla renal", 
                              "2" ~ "Muerto sin falla renal", 
                              NA ~ "Dato perdido"), 
         eventdf = factor(eventdf, 
                          levels = c("Vivo y sin falla renal", 
                                     "Falla renal", 
                                     "Muerto sin falla renal", 
                                     "Total",
                                     "Dato perdido")), 
         time_tipo = case_when(time > 0 & time < 15 ~ "Valor Plausible", 
                               time <= 0 ~ "Valor Implausible", 
                               is.na(time) ~ "Dato perdido", 
                               TRUE ~ as.character(NA))) |> 
  tabyl(time_tipo) |> 
  adorn_pct_formatting() -> tabla_tiempo_perdidos

tabla_tiempo_perdidos |> 
  kbl()

time_tipo	n	percent
Dato perdido	5433	4.0%
Valor Implausible	523	0.4%
Valor Plausible	131519	95.7%

Se aprecia que el porcentaje de datos perdidos de tiempo es de 4.0% (n = 5433). Asimismo, los valores de tiempo implausible son ínfimos y representan el 0.4% (n = 523).

1.5 Seleccion de individuos para el análisis

El total de individuos elegibles es el siguiente:

data_eleg <- data |> 
  filter(age >= 18, ckd_stage == "Stages 3-4")

nrow(data_eleg )

[1] 29293

El número de individuos menores de 18 años es el siguiente:

data_noeleg_age_menos18 <- data |> 
  filter(age < 18)

nrow(data_noeleg_age_menos18)

[1] 136

El número de individuos con datos perdidos en edad es el siguiente:

data_noeleg_age_na <- data |> 
  filter(is.na(age))

nrow(data_noeleg_age_na)

[1] 5468

El número de individuos con diagnostico de CKD diferete a estadio 3a-3:

data_noeleg_ckd12 <- data |> 
  filter(ckd_stage == "Stages 1-2 y 5", eGFR_ckdepi >= 60)

nrow(data_noeleg_ckd12)

[1] 97959

data_noeleg_ckd5 <- data |> 
  filter(ckd_stage == "Stages 1-2 y 5", eGFR_ckdepi < 15)

nrow(data_noeleg_ckd5)

[1] 717

El número de individuos con datos perdidos en el diagnostico de CKD :

data_noeleg_ckd_na <- data |> 
  filter(is.na(ckd_stage))

nrow(data_noeleg_ckd_na)

[1] 9499

El numero de individuos con datos no elegibles de edad o ckd stages es:

data_noeleg_ckd_age <- data |> 
  filter(ckd_stage == "Stages 1-2 y 5" | age < 18)

nrow(data_noeleg_ckd_age)

[1] 98683

El numero de individuos con datos perdidos de edad o perdidos en CKD stages:

data_noeleg_ckd_age_na <- data |> 
  filter(is.na(ckd_stage) | is.na(age))

nrow(data_noeleg_ckd_age_na)

[1] 9499

El numero de individuos potencialmente elgibles por tener edad o CKD stages en el rango o tener datos perdidos

data_eleg_potent <- data |> 
  filter((age >= 18 | is.na(age)) & (ckd_stage == "Stages 3-4" | is.na(ckd_stage)))

nrow(data_eleg_potent)

[1] 38792

El numero individuos incluidos 3a-4:

data_includ <- data_eleg |> 
  filter(time > 0, !is.na(eventd))

nrow(data_includ)

[1] 28843

El numero individuos incluidos 3b-4:

data_includ2 <- data_includ |> 
  filter(ckd_stage2 == "Stages 3b-4")

nrow(data_includ2)

[1] 10596

Los datos excluidos por datos implausibles o perdidos por tiempo son los siguientes

datos_exclud_time_implau <- data_eleg |> 
  filter(time <= 0)

nrow(datos_exclud_time_implau)

[1] 358

Los datos excluidos por datos perdidos en tiempo son los siguientes:

datos_exclud_time_na <- data_eleg |> 
  filter(is.na(time))

nrow(datos_exclud_time_na)

[1] 92

Los datos excluidos por datos perdidos en el status del desenlace:

datos_exclud_eventd_na <- data_eleg |> 
  filter(is.na(eventd))

nrow(datos_exclud_eventd_na)

[1] 92

Los datos excluidos por datos perdidos en el status del desenlace/tiempo o tiempo implausible:

datos_exclud_time_eventd <- data_eleg |> 
  filter(is.na(eventd) | is.na(time) | time <= 0)

nrow(datos_exclud_time_eventd)

[1] 450

1.6 Flujograma de selección / inclusión de participantes

# Create grid of 100 x 100----
data_flow <- tibble(x = 1:100, y = 1:100)

data_flow  %>% 
  ggplot(aes(x, y)) + 
  scale_x_continuous(minor_breaks = seq(1, 100, 1)) + 
  scale_y_continuous(minor_breaks = seq(1, 100, 1)) + 
  theme_linedraw() -> p

# Create boxes----
# 

box_xmin <- 33 - 20
box_xmax <- 75 - 20
box_ymin <- 94
box_ymax <- 99
box_size <- 0.25

text_param <- function(box_min, box_max) {
  mean(c(box_min, box_max))
}

text_size <- 2


p + 
  # Col 0----
  ## Level 1----
  geom_rect(xmin = box_xmin, xmax = box_xmax, 
            ymin = box_ymin - 1, ymax = box_ymax + 1, 
            color = "black", fill = "white", 
            size = box_size) + 
  annotate('text', 
           x = text_param(box_xmin, box_xmax), 
           y = text_param(box_ymin - 1, box_ymax + 1), 
           label = str_glue('Total patients in VISARE database\n(n = {nrow(data)})'), 
           size = text_size ) + 
  ## Level 1.5----
  geom_rect(xmin = box_xmin, xmax = box_xmax, 
            ymin = box_ymin - 27 - 2, ymax = box_ymax - 27 + 2, 
            color = "black", fill = "white", 
            size = box_size) + 
  annotate('text', 
           x = text_param(box_xmin, box_xmax), 
           y = text_param(box_ymin - 27, box_ymax - 27), 
           label = str_glue('Total patients potentially elegible\nwith CKD G3a-G4 and ≥ 18 years old\n(n = {nrow(data_eleg_potent)})'), 
           size = text_size ) + 
  ## Level 2----
  geom_rect(xmin = box_xmin, xmax = box_xmax, 
            ymin = box_ymin - 50 - 1, ymax = box_ymax - 50 + 1, 
            color = "black", fill = "white", 
            size = box_size) + 
  annotate('text', 
           x = text_param(box_xmin, box_xmax), 
           y = text_param(box_ymin - 50 - 1, box_ymax - 50 + 1), 
           label = str_glue('Total patients included in the study\n(n = {nrow(data_eleg)})'), 
           size = text_size ) +   
  # Col -1----
  geom_rect(xmin = box_xmin - 13, xmax = box_xmax - 27, 
            ymin = box_ymin - 81, ymax = box_ymax - 79, 
            color = "black", fill = "white", 
            size = box_size) + 
  annotate('text', 
           x = text_param(box_xmin - 13, box_xmax - 27), 
           y = text_param(box_ymin - 81, box_ymax - 79), 
           label = str_glue('Patients with CKD G3a-G4\n(n = {nrow(data_includ)})'), 
           size = text_size ) + 
  # Col +1----
  ## Level 1----
  geom_rect(xmin = box_xmin + 23, xmax = box_xmax + 47, 
            ymin = box_ymin - 19 + 2.5, ymax = box_ymax - 12 + 2.5, 
            color = "black", fill = "white", 
            size = box_size) + 
  annotate('text', 
           x = text_param(box_xmin + 23, box_xmax + 47), 
           y = text_param(box_ymin - 19 + 2.5, box_ymax - 12 + 2.5), 
           label = str_glue(paste0('Not elegible by age < 18 or CKD in other stages (n = {nrow(data_noeleg_ckd_age)})\n',
                                   '{nrow(data_noeleg_age_menos18)} were less than 18 years old\n', 
                                   '{nrow(data_noeleg_ckd12)} had normal or mildly reduction of eGFR (CKD Stage G1 or G2) \n', 
                                   '{nrow(data_noeleg_ckd5)} had kidney failure (CKD stages G5)\n')), 
           size = text_size )  + 
  ## Level 1.5----
  geom_rect(xmin = box_xmin + 23, xmax = box_xmax + 47, 
            ymin = box_ymin - 19 - 25 + 1.5, ymax = box_ymax - 12 - 25 + 1.5, 
            color = "black", fill = "white", 
            size = box_size) + 
  annotate('text', 
           x = text_param(box_xmin + 23, box_xmax + 47), 
           y = text_param(box_ymin - 19 - 25 + 1.5, box_ymax - 12 - 25 + 1.5), 
           label = str_glue(paste0('Excluded by missing data in selection criteria (n = {nrow(data_noeleg_ckd_age_na)} [{round(100 * nrow(data_noeleg_ckd_age_na)/nrow(data_eleg), 1)}%])\n',
                                   '{nrow(data_noeleg_age_na)} ({round(100 * nrow(data_noeleg_age_na)/nrow(data_eleg), 1)}%) lack of data in age\n', 
                                   '{nrow(data_noeleg_ckd_na)} ({round(100 * nrow(data_noeleg_ckd_na)/nrow(data_eleg), 1)}%) lack of data in eGFR to define CKD stages\n')), 
           size = text_size )  + 
  ## Level 2----
  geom_rect(xmin = box_xmin + 23, xmax = box_xmax + 47, 
            ymin = box_ymin - 19 - 45 - 2, ymax = box_ymax - 12 - 45 - 2, 
            color = "black", fill = "white", 
            size = box_size) + 
  annotate('text', 
           x = text_param(box_xmin + 23, box_xmax + 47), 
           y = text_param(box_ymin - 19 - 45 - 2, box_ymax - 12 - 45 - 2), 
           label = str_glue(paste0('Excluded by missing/implausible data in outcome (n = {nrow(datos_exclud_time_eventd)} [{round(100 * nrow(datos_exclud_time_eventd)/nrow(data_eleg), 1)}%])\n', 
                                   '{nrow(datos_exclud_eventd_na)} ({round(100 * nrow(datos_exclud_eventd_na)/nrow(data_eleg), 1)}%) lack of data in outcome status and time\n', 
                                   '{nrow(datos_exclud_time_implau)} ({round(100 * nrow(datos_exclud_time_implau)/nrow(data_eleg), 1)}%) had negative or zero time-to-event values')), 
           size = text_size )  + 
  ## Level 3----
  geom_rect(xmin = box_xmin + 27, xmax = box_xmax + 13, 
            ymin = box_ymin - 81 , ymax = box_ymax - 79, 
            color = "black", fill = "white", 
            size = box_size) + 
  annotate('text', 
           x = text_param(box_xmin + 27, box_xmax + 13), 
           y = text_param(box_ymin - 81, box_ymax - 79), 
           label = str_glue('Patients with CKD G3b-G4\n(n = {nrow(data_includ2)})'), 
           size = text_size )  + 
  # vertical arrow
  geom_segment(x = text_param(box_xmin, box_xmax), xend = text_param(box_xmin, box_xmax), 
               y = 93, yend = 74, 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt",
               arrow = arrow(length = unit(1, "mm"), type = "closed")) +  
    # vertical arrow
  geom_segment(x = text_param(box_xmin, box_xmax), xend = text_param(box_xmin, box_xmax), 
               y = 65, yend = 50, 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt",
               arrow = arrow(length = unit(1, "mm"), type = "closed")) +  
  geom_segment(x = text_param(box_xmin, box_xmax), xend = text_param(box_xmin, box_xmax), 
               y =43, yend = 25, 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt") + 
  # horizontal arrow 1-->
  geom_segment(x = text_param(box_xmin, box_xmax), xend = box_xmin + 23, 
               y = text_param(box_ymin - 19 + 2, box_ymax - 12 + 2), yend = text_param(box_ymin - 19 + 2, box_ymax - 12 + 2), 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt", 
               arrow = arrow(length = unit(1, "mm"), type = "closed")) + 
  # horizontal arrow 2-->
  geom_segment(x = text_param(box_xmin, box_xmax), xend = box_xmin + 23, 
               y = text_param(box_ymin - 19 - 45, box_ymax - 12 - 45), yend = text_param(box_ymin - 19 - 45, box_ymax - 12 - 45), 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt", 
               arrow = arrow(length = unit(1, "mm"), type = "closed")) + 
  # horizontal arrow 2-->
  geom_segment(x = text_param(box_xmin, box_xmax), xend = box_xmin + 23, 
               y = text_param(box_ymin - 19 - 25 + 1, box_ymax - 12 - 25 + 1), yend = text_param(box_ymin - 19 - 25 + 1, box_ymax - 12 - 25 + 1), 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt", 
               arrow = arrow(length = unit(1, "mm"), type = "closed")) + 
  # horizontal segment --
  geom_segment(x = text_param(box_xmin - 13, box_xmax - 27), xend = text_param(box_xmin + 27, box_xmax + 13), 
               y = 25, yend = 25, 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt") + 
  # vertical arrow -->
  geom_segment(x = text_param(box_xmin - 13, box_xmax - 27), xend = text_param(box_xmin - 13, box_xmax - 27), 
               y = 25, yend = box_ymax - 79, 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt", 
               arrow = arrow(length = unit(1, "mm"), type = "closed")) + 
  # vertical arrow -->
  geom_segment(x = text_param(box_xmin + 27, box_xmax + 13), xend = text_param(box_xmin + 27, box_xmax + 13), 
               y = 25, yend = box_ymax - 79, 
               size = 0.15, 
               linejoin = "mitre", 
               lineend = "butt", 
               arrow = arrow(length = unit(1, "mm"), type = "closed")) + 
  theme_void() +
  theme(plot.background = element_rect(fill = "white")) -> plot_flowchart

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

ggsave(filename = "plot_flowchart.png", 
       plot = plot_flowchart, 
       device = "png", 
       path = here("Figures"), 
       scale = 1, 
       width = 12, 
       height = 12, 
       units = "cm", 
       dpi = 600)

knitr::include_graphics(here("Figures", "plot_flowchart.png"))

1.7 Distribucion de datos según region

data |> 
  mutate(inclusion = case_when(age >= 18 & time > 0 & ckd_stage == "Stages 3-4" ~ "Incluido", 
                               TRUE ~ "No incluido"), 
         cas = fct_rev(fct_infreq(cas))) |> 
  ggplot(aes(x = cas, fill = factor(inclusion, level = c("Incluido", "No incluido")))) + 
  geom_bar() +  
  labs(fill = "Inclusión", x = NULL, y = "Frecuencia absoluta") + 
  coord_flip(expand = FALSE) + 
  scale_y_continuous(limits = c(0, 30000)) + 
  theme_bw() -> p1 

data |> 
  mutate(inclusion = case_when(age >= 18 & time > 0 & ckd_stage == "Stages 3-4" ~ "Incluido", 
                               TRUE ~ "No incluido"), 
         cas = fct_rev(fct_infreq(cas))) |> 
  ggplot(aes(x = cas, fill = factor(inclusion, level = c("Incluido", "No incluido")))) + 
  geom_bar(position = "fill") +  
  labs(fill = "Inclusión", x = NULL, y = "Proporción") + 
  coord_flip(expand = FALSE) + 
  scale_y_continuous(limits = c(0, 1)) + 
  theme_bw()  + 
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) -> p2

data |> 
  mutate(inclusion = case_when(age >= 18 & time > 0 & ckd_stage == "Stages 3-4" ~ "Incluido", 
                               TRUE ~ "No incluido"), 
         cas = fct_rev(fct_infreq(cas))) |> 
  filter(inclusion == "Incluido") |> 
  ggplot(aes(x = cas)) + 
  geom_bar(fill = "#F8766D") +  
  labs(x = NULL, y = "Frecuencia absoluta") + 
  coord_flip(expand = FALSE) + 
  scale_y_continuous(limits = c(0, 11000)) + 
  theme_bw() + 
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) -> p3

(p1 | p3 | p2) + 
  plot_layout(guides = "collect") -> plot_region

ggsave(filename = "plot_region.png", 
       plot = plot_region, 
       device = "png", 
       path = here("Figures"), 
       scale = 2, 
       width = 12,
       height = 9, 
       units = "cm", 
       dpi = 600, 
       bg = "white") 

1.8 Selección de pacientes incluidos en estudio

# Subset patients with CKD Stages 3a-3b-4
data %>%  
  # Criterio de seleccion
  filter(age >= 18, time > 0) |> 
  filter(ckd_stage == "Stages 3-4") |> 
  # Preparacion de datos
  select(cas, sex, age, hta, dm, crea, ckd_stage2, 
         eGFR_ckdepi, acr, urine_album, 
         urine_crea, time5y, eventd5y, 
         grf_cat, acr_cat, ckd_class,
         death2y, eventd2ylab, death5y, 
         eventd5ylab, eventd, time) |> 
  mutate(hta = factor(hta), 
         dm = factor(dm),
         cas2 = case_when(cas %in% c("Lima - Rebagliati") ~ "Lima - Rebagliati", 
                          cas %in% c("Lima - Almenara", "Lima - Sabogal") ~ "Lima Otros", 
                          TRUE ~ "Otras Redes"), 
         cas = fct_rev(fct_infreq(cas)), 
         cas2 = fct_rev(fct_infreq(cas2))
         ) |> 
  mutate(acr = Winsorize(acr, probs = c(0.01, 0.99), na.rm = TRUE), 
         urine_album = Winsorize(urine_album, probs = c(0.01, 0.99), na.rm = TRUE), 
         urine_crea = Winsorize(urine_crea, probs = c(0.01, 0.99), na.rm = TRUE), 
         eventd2ylab = factor(eventd2ylab, 
                              levels = c("Alive w/o Kidney Failure", 
                                         "Kidney Failure", 
                                         "Death w/o Kidney Failure")), 
         eventd5ylab = factor(eventd5ylab, 
                              levels = c("Alive w/o Kidney Failure", 
                                         "Kidney Failure", 
                                         "Death w/o Kidney Failure"))) -> dataA

# Subset patients with CKD Stages 3b-4
data %>% 
  # Criterio de seleccion
  filter(age >= 18, time > 0) |> 
  filter(ckd_stage2 == "Stages 3b-4") |> 
  # Preparacion de datos
  select(cas, sex, age, hta, dm, crea, 
         eGFR_ckdepi, acr, urine_album, 
         urine_crea, time5y, eventd5y, 
         grf_cat, acr_cat, ckd_class,  
         death2y, eventd2ylab, death5y, 
         eventd5ylab, eventd, time) |> 
  mutate(hta = factor(hta), 
         dm = factor(dm),
         cas2 = case_when(cas %in% c("Lima - Rebagliati") ~ "Lima - Rebagliati", 
                          cas %in% c("Lima - Almenara", "Lima - Sabogal") ~ "Lima Otros", 
                          TRUE ~ "Otras Redes"), 
         cas = fct_rev(fct_infreq(cas)), 
         cas2 = fct_rev(fct_infreq(cas2))
         ) |> 
  mutate(acr = Winsorize(acr, probs = c(0.01, 0.99), na.rm = TRUE), 
         urine_album = Winsorize(urine_album, probs = c(0.01, 0.99), na.rm = TRUE), 
         urine_crea = Winsorize(urine_crea, probs = c(0.01, 0.99), na.rm = TRUE), 
         eventd2ylab = factor(eventd2ylab, 
                              levels = c("Alive w/o Kidney Failure", 
                                         "Kidney Failure", 
                                         "Death w/o Kidney Failure")), 
         eventd5ylab = factor(eventd5ylab, 
                              levels = c("Alive w/o Kidney Failure", 
                                         "Kidney Failure", 
                                         "Death w/o Kidney Failure"))) -> dataB

1.9 Cargar Funciones Escritas por Usuario

source(here("Code", "source", "kfre_pi.R"))
source(here("Code", "source", "kfre_pr.R"))
source(here("Code", "source", "oe_ratio.R"))
source(here("Code", "source", "calibration_intercept.R"))
source(here("Code", "source", "calibration_slope.R"))
source(here("Code", "source", "auc.R"))
source(here("Code", "source", "validate.mids.R"))
source(here("Code", "source", "pool.validate.mids.R"))
source(here("Code", "source", "pool.auc.mids.R"))
source(here("Code", "source", "predict.mira.R"))
source(here("Code", "source", "process_imp_cal_plot.R"))

1.10 Creacion de hazard acumulado causa especifico mediante estimador de Nelson-Aalen

cumhaz1 <- basehaz(coxph(Surv(time, eventd == 1) ~ 1, data = dataA)) |> 
  rename(cumhaz1 = hazard)
cumhaz2 <- basehaz(coxph(Surv(time, eventd == 2) ~ 1, data = dataA))|> 
  rename(cumhaz2 = hazard)

dataA <- dataA |> 
  left_join(cumhaz1, by = "time") |> 
  left_join(cumhaz2, by = "time")

1.11 Preparación de datos

dataA_imp <- dataA |> 
  mutate(eventdf = factor(eventd, levels = c(0, 1, 2))) |> 
  select(-ckd_class, -acr_cat) |> 
  mutate(sex_cumhaz1 = as.integer(sex == "Masculino") * (cumhaz1 - mean(cumhaz1)), 
         age_cumhaz1 = (age - mean(age)) * (cumhaz1 - mean(cumhaz1)), 
         eGFR_ckdepi_cumhaz1 = (eGFR_ckdepi - mean(age)) * (cumhaz1 - mean(cumhaz1)), 
         sex_cumhaz2 = as.integer(sex == "Masculino") * (cumhaz2 - mean(cumhaz2)), 
         age_cumhaz2 = (age - mean(age)) * (cumhaz2 - mean(cumhaz2)), 
         eGFR_ckdepi_cumhaz2 = (eGFR_ckdepi - mean(age)) * (cumhaz2 - mean(cumhaz2)), 
         log_urine_crea = log(urine_crea), 
         log_urine_album = log(urine_album),
         log_acr = log(acr)) |> 
  select(-urine_crea, -urine_album, -acr) |> 
  relocate(all_of(c("hta", "dm", "log_urine_crea", "log_acr", "log_urine_album")), 
           .after = eGFR_ckdepi_cumhaz2)

1.12 Exploración de datos perdidos

Las variables a tener en cuenta se muestran a continuacion:

dataA |> 
  glimpse()

Rows: 28,843
Columns: 25
$ cas         <fct> Piura, Piura, Cusco, Cusco, Cusco, Cusco, Cusco, Cusco, Cu…
$ sex         <fct> Femenino, Masculino, Femenino, Femenino, Masculino, Femeni…
$ age         <dbl> 99, 100, 99, 99, 99, 99, 99, 99, 99, 100, 99, 99, 99, 96, …
$ hta         <fct> 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0…
$ dm          <fct> 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1…
$ crea        <dbl> 1.00, 1.03, 0.96, 1.13, 1.31, 0.91, 0.98, 1.37, 1.30, 1.15…
$ ckd_stage2  <fct> Stages 1-3 y 5, Stages 1-3 y 5, Stages 1-3 y 5, Stages 3b-…
$ eGFR_ckdepi <dbl> 46.67256, 59.33391, 49.03382, 40.26148, 44.67792, 52.30950…
$ acr         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ urine_album <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ urine_crea  <dbl> NA, 120.87, NA, 62.80, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ time5y      <dbl> 5.000000, 5.000000, 5.000000, 5.000000, 5.000000, 5.000000…
$ eventd5y    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 2, 0, 2…
$ grf_cat     <fct> G3a, G3a, G3a, G3b, G3b, G3a, G3a, G3b, G3a, G3a, G3a, G3a…
$ acr_cat     <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ ckd_class   <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ death2y     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0…
$ eventd2ylab <fct> Alive w/o Kidney Failure, Alive w/o Kidney Failure, Alive …
$ death5y     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1…
$ eventd5ylab <fct> Alive w/o Kidney Failure, Alive w/o Kidney Failure, Alive …
$ eventd      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 2, 0, 2…
$ time        <dbl> 9.993155, 9.563313, 9.982204, 8.974675, 8.974675, 8.974675…
$ cas2        <fct> Otras Redes, Otras Redes, Otras Redes, Otras Redes, Otras …
$ cumhaz1     <dbl> 0.06149899, 0.06149899, 0.06149899, 0.06034024, 0.06034024…
$ cumhaz2     <dbl> 0.54808445, 0.51587330, 0.54808445, 0.47474466, 0.47474466…

La cantidad de datos perdidos por variable es la siguiente:

Perdidos por variable

dataA |> 
  select(-c(ckd_stage2, time5y, eventd5y, acr_cat, 
            death2y, eventd2ylab, death5y, eventd5ylab, time, cas)) |> 
  miss_var_summary() |> 
  kbl()

variable	n_miss	pct_miss
urine_album	21358	74.04916
acr	20299	70.37756
ckd_class	20299	70.37756
urine_crea	16851	58.42319
dm	7302	25.31637
hta	4376	15.17179
sex	0	0.00000
age	0	0.00000
crea	0	0.00000
eGFR_ckdepi	0	0.00000
grf_cat	0	0.00000
eventd	0	0.00000
cas2	0	0.00000
cumhaz1	0	0.00000
cumhaz2	0	0.00000

Patrón

A continuación se muestra un gráfico de patrón de datos perdidos. Podemos apreciar un patrón de datos general, caracterizado por tener datos perdidos multivariado, no se aprecia ningún patrón monótono de datos perdidos (como era de esperarse) y tampoco existe file matching. ASimismo, se aprecia que los datos perdidos muestran un patrón conectado a nivel de filas y columnas.

dataA |> 
  select(-c(ckd_stage2, time5y, eventd5y, acr_cat, ckd_class,
            death2y, eventd2ylab, death5y, eventd5ylab, time, cas)) |> 
  plot_pattern(rotate = TRUE) -> plot_patterns

ggsave(filename = "plot_patterns.png", 
       plot = plot_patterns, 
       device = "png", 
       path = here("Figures", "Imputation_Diagnostic"), 
       scale = 2, 
       width = 9,
       height = 9, 
       units = "cm", 
       dpi = 600, 
       bg = "white") 

1.13 Influx - outflux

dataA |> 
  select(-c(ckd_stage2, time5y, eventd5y, acr_cat, ckd_class,
            death2y, eventd2ylab, death5y, eventd5ylab, cas)) |> 
  plot_flux(label = FALSE) -> plot_influx

# Primero, averiguamos cuántas capas hay.
num_layers <- length(plot_influx$layers)

# Examinamos cada capa para encontrar la geom_point() que no deseamos.
# Esto imprimirá las capas y deberías buscar la que contiene la geom_point sin shape.
for (i in 1:num_layers) {
  print(plot_influx$layers[[i]])
}

mapping: intercept = ~intercept, slope = ~slope 
geom_abline: na.rm = FALSE
stat_identity: na.rm = FALSE
position_identity 
geom_point: na.rm = FALSE
stat_identity: na.rm = FALSE
position_nudge 

# Eliminamos la capa geom_point que no queremos.
# La salida muestra que es la segunda capa, así que la eliminamos.
plot_influx$layers <- plot_influx$layers[-2] 

# Asegúrate de tener suficientes formas para cada nivel único de la variable.
unique_vrbs <- unique(plot_influx$data$vrb)
shapes <- seq_along(unique_vrbs)

# Ahora deberías volver a agregar la capa geom_point con las formas y colores adecuados.
plot_influx <- plot_influx + 
  geom_jitter(aes(shape = vrb, colour = vrb), width = 0.025, height = 0.025) +
  scale_shape_manual(values = shapes) +
  guides(colour = guide_legend(override.aes = list(shape = shapes)),
         shape = FALSE)

Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
of ggplot2 3.3.4.

ggsave(filename = "plot_influx.png", 
       plot = plot_influx, 
       device = "png", 
       path = here("Figures", "Imputation_Diagnostic"), 
       scale = 2, 
       width = 9,
       height = 9, 
       units = "cm", 
       dpi = 600, 
       bg = "white") 

1.14 Por region

dataA |> 
  mutate(miss_acr = if_else(is.na(acr), "Missing", "Completo"), 
         miss_acr = factor(miss_acr, level = c("Completo", "Missing"))) |> 
  ggplot(aes(x = cas, fill = miss_acr)) + 
  geom_bar(position = "fill") +  
  labs(fill = "Datos perdidos", x = NULL, y = "Proporción") + 
  coord_flip(expand = FALSE) + 
  scale_y_continuous(limits = c(0, 1)) + 
  theme_bw() -> p1

dataA |> 
  mutate(miss_acr = if_else(is.na(acr), "Missing", "Completo"), 
         miss_acr = factor(miss_acr, level = c("Completo", "Missing"))) |> 
  ggplot(aes(x = cas, fill = miss_acr)) + 
  geom_bar() +  
  labs(fill = "Datos perdidos", x = NULL, y = "Frecuencia absoluta") + 
  coord_flip(expand = FALSE) + 
  scale_y_continuous(limits = c(0, 10000)) + 
  theme_bw()  -> p2

(p1 | p2) + plot_layout(guides = "collect") -> plot_region_missing

ggsave(filename = "plot_region_missing.png", 
       plot = plot_region_missing, 
       device = "png", 
       path = here("Figures", "Imputation_Diagnostic"), 
       scale = 2, 
       width = 16,
       height = 8, 
       units = "cm", 
       dpi = 600, 
       bg = "white") 

dataA |> 
  mutate(miss_acr = if_else(is.na(acr), "Missing", "Completo"), 
         miss_acr = factor(miss_acr, level = c("Completo", "Missing"))) |> 
  ggplot(aes(x = cas2, fill = miss_acr)) + 
  geom_bar(position = "fill") +  
  labs(fill = "Datos perdidos", x = NULL, y = "Proporción") + 
  coord_flip(expand = FALSE) + 
  scale_y_continuous(limits = c(0, 1)) + 
  theme_bw() -> p1

dataA |> 
  mutate(miss_acr = if_else(is.na(acr), "Missing", "Completo"), 
         miss_acr = factor(miss_acr, level = c("Completo", "Missing"))) |> 
  ggplot(aes(x = cas2, fill = miss_acr)) + 
  geom_bar() +  
  labs(fill = "Datos perdidos", x = NULL, y = "Frecuencia absoluta") + 
  coord_flip(expand = FALSE) + 
  scale_y_continuous(limits = c(0, 16000)) + 
  theme_bw()  -> p2

(p1 | p2) + plot_layout(guides = "collect") -> plot_region_missing

ggsave(filename = "plot_region_missing2.png", 
       plot = plot_region_missing, 
       device = "png", 
       path = here("Figures", "Imputation_Diagnostic"), 
       scale = 2, 
       width = 16,
       height = 6, 
       units = "cm", 
       dpi = 600, 
       bg = "white") 

1.15 Inicializar MICE

Se inicia con maxit = 0 para verificar que no haya errores

setup_impA <- mice(dataA_imp, maxit = 0)

Warning: Number of logged events: 2

Se verifica si hay eventos que revisar

setup_impA$loggedEvents

Se revisa que modelos de imputacion

make.method(dataA_imp)

                cas                 sex                 age                crea 
                 ""                  ""                  ""                  "" 
         ckd_stage2         eGFR_ckdepi              time5y            eventd5y 
                 ""                  ""                  ""                  "" 
            grf_cat             death2y         eventd2ylab             death5y 
                 ""                  ""                  ""                  "" 
        eventd5ylab              eventd                time                cas2 
                 ""                  ""                  ""                  "" 
            cumhaz1             cumhaz2             eventdf         sex_cumhaz1 
                 ""                  ""                  ""                  "" 
        age_cumhaz1 eGFR_ckdepi_cumhaz1         sex_cumhaz2         age_cumhaz2 
                 ""                  ""                  ""                  "" 
eGFR_ckdepi_cumhaz2                 hta                  dm      log_urine_crea 
                 ""            "logreg"            "logreg"               "pmm" 
            log_acr     log_urine_album 
              "pmm"               "pmm" 

predA <- make.predictorMatrix(dataA_imp)
plot_pred(predA, rotate = TRUE) 

predA[, "crea"] <- 0
predA[, "time5y"] <- 0
predA[, "eventd5y"] <- 0
predA[, "death2y"] <- 0
predA[, "eventd2ylab"] <- 0
predA[, "death5y"] <- 0
predA[, "eventd5ylab"] <- 0
predA[, "eventd"] <- 0
predA[, "time"] <- 0
predA[, "ckd_stage2"] <- 0
predA[, "grf_cat"] <- 0
predA[, "cas"] <- 0


predA["crea", ] <- 0
predA["time5y", ] <- 0
predA["eventd5y", ] <- 0
predA["death2y", ] <- 0
predA["eventd2ylab", ] <- 0
predA["death5y", ] <- 0
predA["eventd5ylab", ] <- 0
predA["eventd", ] <- 0
predA["time", ] <- 0
predA["ckd_stage2", ] <- 0
predA["grf_cat", ] <- 0
predA["cas", ] <- 0

predA[c("log_urine_crea", "log_acr"), c("log_urine_crea", "log_acr")] <- 0
predA[c("log_urine_album", "log_acr"), c("log_urine_album", "log_acr")] <- 0

methA <- make.method(dataA_imp)
methA

                cas                 sex                 age                crea 
                 ""                  ""                  ""                  "" 
         ckd_stage2         eGFR_ckdepi              time5y            eventd5y 
                 ""                  ""                  ""                  "" 
            grf_cat             death2y         eventd2ylab             death5y 
                 ""                  ""                  ""                  "" 
        eventd5ylab              eventd                time                cas2 
                 ""                  ""                  ""                  "" 
            cumhaz1             cumhaz2             eventdf         sex_cumhaz1 
                 ""                  ""                  ""                  "" 
        age_cumhaz1 eGFR_ckdepi_cumhaz1         sex_cumhaz2         age_cumhaz2 
                 ""                  ""                  ""                  "" 
eGFR_ckdepi_cumhaz2                 hta                  dm      log_urine_crea 
                 ""            "logreg"            "logreg"               "pmm" 
            log_acr     log_urine_album 
              "pmm"               "pmm" 

plot_pred(predA,  
          rotate = TRUE) -> plot_matriz_pred

ggsave(filename = "plot_matriz_pred.png", 
       plot = plot_matriz_pred, 
       device = "png", 
       path = here("Figures", "Imputation_Diagnostic"), 
       scale = 2, 
       width = 9,
       height = 9, 
       units = "cm", 
       dpi = 600, 
       bg = "white") 

1.16 Imputation

n_cores <- detectCores()
n_cores

[1] 8

# Evaluate futures in parallel - max of two workers to avoid hogging resources
future::plan("multisession", workers = n_cores)

tic()
data_impA <- futuremice(dataA_imp, 
                        method = methA, 
                        m = 100, 
                        maxit = 500, 
                        pred = predA, 
                        n.core = n_cores)
toc()

13421.53 sec elapsed

data_impA$loggedEvents

NULL

semilla_reproducible <- data.frame(semilla = data_impA$parallelseed)

data_impA$visitSequence

 [1] "cas"                 "sex"                 "age"                
 [4] "crea"                "ckd_stage2"          "eGFR_ckdepi"        
 [7] "time5y"              "eventd5y"            "grf_cat"            
[10] "death2y"             "eventd2ylab"         "death5y"            
[13] "eventd5ylab"         "eventd"              "time"               
[16] "cas2"                "cumhaz1"             "cumhaz2"            
[19] "eventdf"             "sex_cumhaz1"         "age_cumhaz1"        
[22] "eGFR_ckdepi_cumhaz1" "sex_cumhaz2"         "age_cumhaz2"        
[25] "eGFR_ckdepi_cumhaz2" "hta"                 "dm"                 
[28] "log_urine_crea"      "log_acr"             "log_urine_album"    

plot_trace(data_impA) + 
  theme(legend.position="none") -> plot_trace_dx

ggsave(filename = "plot_trace_dx.png", 
       plot = plot_trace_dx, 
       device = "png", 
       path = here("Figures", "Imputation_Diagnostic"), 
       scale = 2.5, 
       width = 12,
       height = 6, 
       units = "cm", 
       dpi = 600, 
       bg = "white") 

mean_converg <- convergence(data_impA, diagnostic = "all", parameter = "mean") |> 
  mutate(type = "mean")
sd_converg <- convergence(data_impA, diagnostic = "all", parameter = "sd") |> 
  mutate(type = "sd") 

mean_converg |> 
  bind_rows(sd_converg) |> 
  filter(vrb %in% c("hta", "dm", "log_urine_crea", "log_urine_album", 
                    "log_acr")) |> 
  ggplot(aes(x = .it, y = ac, color = type)) + 
  geom_line() + 
  geom_hline(yintercept = 0) + 
  facet_grid(type ~ vrb) + 
  theme_bw() -> plot_autocor_dx

ggsave(filename = "plot_autocor_dx.png", 
       plot = plot_autocor_dx, 
       device = "png", 
       path = here("Figures", "Imputation_Diagnostic"), 
       scale = 2.5, 
       width = 12,
       height = 6, 
       units = "cm", 
       dpi = 600, 
       bg = "white") 

Warning: Removed 2 rows containing missing values (`geom_line()`).

mean_converg |> 
  bind_rows(sd_converg) |> 
  filter(vrb %in% c("hta", "dm", "log_urine_crea", "log_urine_album", 
                    "log_acr")) |> 
  ggplot(aes(x = .it, y = psrf, color = type)) + 
  geom_line() + 
  geom_hline(yintercept = 1) + 
  facet_grid(type ~ vrb) + 
  theme_bw()-> plot_psrf_dx

ggsave(filename = "plot_psrf_dx.png", 
       plot = plot_psrf_dx, 
       device = "png", 
       path = here("Figures", "Imputation_Diagnostic"), 
       scale = 2.5, 
       width = 12,
       height = 6, 
       units = "cm", 
       dpi = 600, 
       bg = "white") 

densityplot(data_impA, ~ log_urine_crea)

Hint: Did you know, an equivalent figure can be created with `ggmice()`?
For example, to plot a variable named 'my_vrb' from a mids object called
'my_mids', run:
  ggmice(my_mids, ggplot2::aes(x = my_vrb, group = .imp)) +
  ggplot2::geom_density()
ℹ See amices.org/ggmice for more info.
This message is displayed once per session.

densityplot(data_impA, ~ log_urine_album)

densityplot(data_impA, ~ log_acr)

1.16.1 Datos apilados

saveRDS(data_impA, here("Data", "Tidy", "data_impA.rds"))
saveRDS(semilla_reproducible, here("Data", "Tidy", "semilla_reproducible.rds"))

1.17 Ticket de Reprocubilidad

sessionInfo()

R version 4.3.3 (2024-02-29 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default


locale:
[1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8   
[3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C                  
[5] LC_TIME=Spanish_Spain.utf8    

time zone: America/Lima
tzcode source: internal

attached base packages:
 [1] parallel  grid      splines   stats     graphics  grDevices utils    
 [8] datasets  methods   base     

other attached packages:
 [1] smplot2_0.1.0              impstep_0.1.0             
 [3] ggmice_0.1.0               tictoc_1.2                
 [5] furrr_0.3.1                future_1.33.2             
 [7] ggExtra_0.10.1             gtools_3.9.5              
 [9] DescTools_0.99.54          naniar_1.0.0              
[11] mice_3.16.0                PerformanceAnalytics_2.0.4
[13] xts_0.13.1                 zoo_1.8-12                
[15] VIM_6.2.2                  colorspace_2.1-0          
[17] janitor_2.2.0              gt_0.10.0                 
[19] DiagrammeR_1.0.10          ggtext_0.1.2              
[21] htmltools_0.5.7            devtools_2.4.5            
[23] usethis_2.2.2              gghalves_0.1.4            
[25] xml2_1.3.5                 downlit_0.4.3             
[27] broom_1.0.5                dcurves_0.4.0             
[29] glue_1.6.2                 labelled_2.12.0           
[31] scales_1.3.0               cowplot_1.1.2             
[33] ggsci_3.0.0                survminer_0.4.9           
[35] ggpubr_0.6.0               patchwork_1.1.3           
[37] webshot_0.5.5              gridExtra_2.3             
[39] rsample_1.2.0              lubridate_1.9.3           
[41] forcats_1.0.0              stringr_1.5.1             
[43] dplyr_1.1.4                purrr_1.0.2               
[45] readr_2.1.5                tidyr_1.3.1               
[47] tibble_3.2.1               ggplot2_3.4.4             
[49] tidyverse_2.0.0            boot_1.3-29               
[51] gtsummary_1.7.2            flextable_0.9.4           
[53] kableExtra_1.3.4           knitr_1.45                
[55] plotrix_3.8-4              pec_2023.04.12            
[57] prodlim_2023.08.28         pseudo_1.4.3              
[59] geepack_1.3.9              KMsurv_0.1-5              
[61] mstate_0.3.2               riskRegression_2023.09.08 
[63] cmprsk_2.2-11              rms_6.7-1                 
[65] Hmisc_5.1-1                survival_3.5-8            
[67] skimr_2.1.5                here_1.0.1                

loaded via a namespace (and not attached):
  [1] matrixStats_1.1.0       fs_1.6.3                httr_1.4.7             
  [4] RColorBrewer_1.1-3      repr_1.1.6              numDeriv_2016.8-1.1    
  [7] profvis_0.3.8           tools_4.3.3             backports_1.4.1        
 [10] utf8_1.2.4              R6_2.5.1                jomo_2.7-6             
 [13] urlchecker_1.0.1        withr_3.0.0             sp_2.0-0               
 [16] quantreg_5.97           cli_3.6.1               textshaping_0.3.7      
 [19] pacman_0.5.1            officer_0.6.3           sandwich_3.0-2         
 [22] labeling_0.4.3          mvtnorm_1.2-4           survMisc_0.5.6         
 [25] robustbase_0.99-0       polspline_1.1.24        proxy_0.4-27           
 [28] askpass_1.2.0           QuickJSR_1.0.8          StanHeaders_2.26.28    
 [31] systemfonts_1.0.5       foreign_0.8-86          gfonts_0.2.0           
 [34] svglite_2.1.2           parallelly_1.37.1       sessioninfo_1.2.2      
 [37] pwr_1.3-0               readxl_1.4.3            rstudioapi_0.15.0      
 [40] httpcode_0.3.0          shape_1.4.6.1           visNetwork_2.1.2       
 [43] generics_0.1.3          car_3.1-2               zip_2.3.0              
 [46] inline_0.3.19           loo_2.6.0               Matrix_1.6-1           
 [49] fansi_1.0.6             abind_1.4-5             lifecycle_1.0.4        
 [52] multcomp_1.4-25         yaml_2.3.7              snakecase_0.11.1       
 [55] carData_3.0-5           promises_1.2.1          mitml_0.4-5            
 [58] crayon_1.5.2            miniUI_0.1.1.1          lattice_0.22-5         
 [61] haven_2.5.4             pillar_1.9.0            gld_2.6.6              
 [64] future.apply_1.11.0     codetools_0.2-19        pan_1.9                
 [67] V8_4.4.0                fontLiberation_0.1.0    data.table_1.14.8      
 [70] broom.helpers_1.14.0    remotes_2.4.2.1         vcd_1.4-11             
 [73] png_0.1-8               vctrs_0.6.4             cellranger_1.1.0       
 [76] gtable_0.3.4            cachem_1.0.8            xfun_0.40              
 [79] mime_0.12               iterators_1.0.14        lava_1.7.3             
 [82] ellipsis_0.3.2          TH.data_1.1-2           nlme_3.1-164           
 [85] fontquiver_0.2.1        rstan_2.32.3            rprojroot_2.0.4        
 [88] rpart_4.1.23            nnet_7.3-19             Exact_3.2              
 [91] tidyselect_1.2.1        sdamr_0.2.0             compiler_4.3.3         
 [94] curl_5.1.0              glmnet_4.1-8            rvest_1.0.3            
 [97] htmlTable_2.4.2         expm_0.999-7            SparseM_1.81           
[100] fontBitstreamVera_0.1.1 checkmate_2.3.0         DEoptimR_1.1-3         
[103] lmtest_0.9-40           quadprog_1.5-8          digest_0.6.33          
[106] minqa_1.2.6             rmarkdown_2.25          pkgconfig_2.0.3        
[109] base64enc_0.1-3         lme4_1.1-35.2           highr_0.10             
[112] fastmap_1.1.1           rlang_1.1.2             htmlwidgets_1.6.3      
[115] shiny_1.8.0             farver_2.1.1            jsonlite_1.8.8         
[118] magrittr_2.0.3          Formula_1.2-5           munsell_0.5.0          
[121] Rcpp_1.0.11             gdtools_0.3.4           visdat_0.6.0           
[124] stringi_1.7.12          rootSolve_1.8.2.4       MASS_7.3-60.0.1        
[127] pkgbuild_1.4.4          listenv_0.9.1           lmom_3.0               
[130] mets_1.3.2              gridtext_0.1.5          hms_1.1.3              
[133] timereg_2.0.5           uuid_1.1-1              ranger_0.15.1          
[136] ggsignif_0.6.4          stats4_4.3.3            pkgload_1.3.4          
[139] crul_1.4.0              evaluate_0.23           RcppParallel_5.1.7     
[142] laeken_0.5.2            nloptr_2.0.3            tzdb_0.4.0             
[145] foreach_1.5.2           httpuv_1.6.12           MatrixModels_0.5-3     
[148] openssl_2.1.1           km.ci_0.5-6             xtable_1.8-4           
[151] e1071_1.7-13            rstatix_0.7.2           later_1.3.1            
[154] viridisLite_0.4.2       class_7.3-22            ragg_1.2.6             
[157] memoise_2.0.1           cluster_2.1.5           timechange_0.2.0       
[160] globals_0.16.3